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Abstract. In this paper we provide a detailed convergence analysis for an unconditionally energy 
stable, second-order accurate convex splitting scheme for the Modified Phase Field Crystal equation, 
a generalized damped wave equation for which the usual Phase Field Crystal equation is a special 
degenerate case. The fully discrete, fully second-order finite difference scheme in question was derived 
in a recent work [2]. An introduction of a new variable if), corresponding to the temporal derivative of 
the phase variable <f>, could bring an accuracy reduction in the formal consistency estimate, because 
of the hyperbolic nature of the equation. A higher order consistency analysis by an asymptotic 
expansion is performed to overcome this difficulty. In turn, second order convergence in both time 
and space is established in a discrete L°° (0, T; H 3 ) norm. 
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1. Introduction. The modified phase field crystal (MPFC) equation is given 

by na 

(3d tt <P + d t 4> = A (<p 3 + a(f> + 2Acf> + A 2 <p) , (1.1) 

where j3 > and a > 0. Equation (jl.lj) is a generalized damped wave equation. The 
parabolic phase field crystal (PFC) equation is recovered in the degenerate case when 
P = 0. See [21 [TH HU HU [20] and references therein for the physical motivation for 
the MPFC equation. The existence and uniqueness of global smooth solutions of the 
MPFC equation was established in our recent article [20], assuming that the initial 
data are smooth. Very recently, we devised and implemented a second-order convex 
splitting scheme for the MPFC equation 2] . The solver for the discrete equations was 
based on a nearly optimally efficient nonlinear multigrid method. While we proved a 
priori unconditional stability and unconditional solvability results for the scheme, we 
did not perform a convergence analysis. The goal of this paper is to provide a detailed 
convergence analysis of the second order convex-splitting scheme for MPFC equation 
(jl.ip proposed in [2]. To our knowledge no second order convergence analysis exists 
for scheme for either the PFC or the MPFC equation. 

Because of the close relationship between the MPFC and PFC models, methods 
for the latter equation can be adapted and applied to the former. See, for example, [T] 
|6] [7] [TO] [T3] [22] for some recent approximation methods specifically for the PFC model. 
Methods specifically designed for the MPFC equation can be found in [2[ [T2l [16j [20l 
119] . Stefanovic et al., [16j employed a semi- implicit finite difference discretization, 
with a multigrid algorithm for solving the algebraic equations. They provide no 
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numerical analysis for their scheme, which is significantly different from schemes we 
propose and analyze. The MPFC scheme in [TJ] is more or less the same as the 
first-order convex-splitting that we devised earlier in j^Ol H] • 

The MPFC equation may be viewed as a perturbed gradient flow with 

respect to an energy. Specifically, consider a dimensionless spatial energy of the 
form [8j QI] 

m = L{^ 4+ ^ 2 _ lw|2 + 5 (A0)2 } dx ' (L2) 

where <\> '■ n C R 2 — >• R is the "atom" density field, and a > is a constant. Suppose 
that f2 = (0, Ljc) x (0, L y ) and </i is periodic on ft. Define /z to be the chemical potential 
with respect to E: 

ft := 6^E = (j> 3 + (1 - e)4> + 2A0 + A 2 0, (1.3) 

where (S^i? denotes the variational derivative with respect to <f>. Clearly, the MPFC 
equation may be redefined as 

Pd tt <l> + d t 4> = A/i, (1.4) 

where f3 > 0. As mentioned, when (5 = the PFC equation is recovered. Herein we 
will restrict ourselves to the case that /3 > to avoid degeneracy. See the discussion 
in [5D] for some equations in the literature that are closely related to (ll.lj) . 

First, note that the energy (|1.2p is not necessarily non-increasing in time along 
the solution trajectories of Eq. (I1.4[) . However, solutions of the MPFC equation do 
dissipate a pseudo energy, as we show momentarily. Also observe that Eq. (|1.4p is not 
precisely a mass conservation equation due to the term f3dtt4>- However, it is easy to 

show that if / 9t0(x, 0) dx = 0, then / 9 t </>dx = for all time [501 US]- Herein we 

Ju Jn 
assume 9 t 0(x, 0) = 0, for simplicity, which trivially satisfies the condition for mass 

conservation. 

We now recast the MPFC equation (| 1 .4[) as the following system of equations: 

pdti/> = Afx — tp, d t <j> = i/>. (1.5) 
And we introduce the pseudo energy 

S^^-E^ + ^ml^. (1.6) 

See [19l [20] for precise definitions of the H~ 1 inner product and norm. For well- 
definedness of the H^ 1 norm, we requires that J n tp cht = 0. This is the case since we 
use the initial data 

</>(•, 0) =d t cj>(-,0) = infi. (1.7) 

A simple calculation [19j [20] shows that sufficiently regular solutions dissipate the 
pseudo-energy at the rate 

dtS = -(VsV)jj-i <0. (1.8) 

In other words, the pseudo energy is non-increasing in time. The primary motivation 
in the convex splitting framework is to design fully and semi-discrete schemes that 
mimic this pseudo-energy dissipation [31 032 [5D] . 



SECOND ORDER SCHEME FOR THE MPFC EQUATION 



3 



The first order convex splitting scheme for (jl.ip was proposed and analyzed in a 
recent article (TH], as we have mentioned. However, the extension to the second order 
convergence analysis is highly non-trivial, mainly due to an 0(s 2 ) numerical error 
between the centered difference of cj> and the mid-point average of ip. As observed in 
[2] , the introduction of the variable ip greatly facilitates the numerical implementation. 
However, if one is not careful, the above-mentioned 0(s 2 ) numerical error might seem 
to introduce a reduction of temporal accuracy, because of the second order time 
derivative involved in the equation. To overcome this difficulty in the paper, we 
have to perform a higher order consistency analysis by an asymptotic expansion; as 
a result, the constructed approximate solution satisfies the numerical scheme with a 
higher order truncation error. A projection of the exact solution onto the Fourier 
space is taken so that an optimal regularity requirement is obtained. 

Second order convergence analysis has always been very challenging for nonlinear 
hyperbolic equation with second order temporal derivative involved. The nonlinear 
error term must be carefully expanded, and a discrete Sobolev inequality is needed to 
bound the discrete L°° and W 1 ' norms of the numerical error function. In addition, 
we need to take inner product with the error equation by the (discrete) time derivative 
of the numerical error, because of the hyperbolic nature of MPFC equation. In the 
end, a full second order convergence in a discrete L°° (0,T;H 3 ) norm is established. 

In Sec. [2] we define the second-order convex splitting scheme and restate some 
solvability and stability results from [2j. In Sec. [3] we present the convergence analysis 
for the second order scheme. We give some concluding remarks in Sec. |U Moreover, 
some technical details of the forthcoming analysis are provided in two appendices. 
In App. [A] we give the finite difference background for the analysis, including our 
notation, some of the necessary difference operators, and the some useful inequalities. 
In a second appendix, App. [Bj we give the details of the consistency analyses related 
to our scheme. 

2. The Second-Order Scheme and its Properties. Here we redefine our 
second-order convex splitting scheme from [2J. We also restate some of the uncon- 
ditional solvability and stability results for this scheme. We note that we used a 
different non-dimensional scaling of the MPFC equation (11.11) in [2] than we do here, 
and some of the restated results below will be in a slightly modified form. However, 
this difference is only superficial. The reader is directed to App. [A] for an introduc- 
tion to the notation, as well as some of the standard tools from cell centered finite 
differences, that are used below. 

2.1. Discrete Energy and the Convex-Splitting Scheme. We first intro- 
duce a fully discrete energy that is consistent with the continuous space energy (|1.2[) . 
In particular, define the discrete energy F : Cmxn — > K to be 

m ■■= \ ut + f ml - \\Vh<t>\\l + \ w^t • (2- 1 ) 

The discrete analogue to (jl.6[) is 

F{4>^):=F^) + ^ |H| 2 _ 1; (2.2) 

defined for any </) G C^xn and any ip G H. The norms above, including the "— 1" 
norm, are defined in App. [A] 
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Note that if <f> € Cmxn is periodic, then it is easy to see that the energies 

F c (ct>) = ~IMl4 + f ||^ + ~||A fc <^ and F e (0) = ||V h 0|£ (2.3) 

are convex [T9], [22] . Hence F, as defined in (|2.1j) . admits the convex splitting F = 
F c — F e . Our second-order scheme will exploit this decomposition of F. Eyre [9] is 
often credited with popularizing the idea that the numerical scheme should respect 
the convexity structure of the energy for the purposes of numerical stability and 
solvability. His original scheme was first-order accurate in time and was restricted to 
non-conserved gradient flows. But this approach has been extended to craft schemes 
for a number of gradient-flow equations of parabolic type; see for example HH [T8l [22] 
[2Tj . The convex splitting framework was extended for the hyperbolic MPFC equation 
in P3JJ [5D] . We extended the framework for second-order schemes in [51 [101 03] ■ 
The following second-order convex splitting scheme for the MPFC equation is from 
our recent paper [2]: given </> fc ~ 1 , <j) k , ip k G Cmxn periodic, find </> fc+1 , ip k+1 , [i k+1 / 2 e 
Cmxn periodic such that 

(3 _ ^fc) = s A hl i k+1/2 - s V fc+1/2 , (2-4) 

//' +1/2 = x (0 fe+1 , /) + a<j> k+l l 2 + 2A h $ k+1 ^ 2 + A 2 ^ 1 / 2 , (2.5) 

<P k+l -<p k = s^ k+l / 2 , (2.6) 



where 



2 ' 2 r ' r 2 

It is obvious that %((/>, <fi) = 3 . In [2], we used the initial data 

r 1 ^ , 1&° = 0. (2.7) 

Note that <fr x = cj) is an 0(s 2 ) approximation to the phase variable at the time 
"ghost" point k = — 1; such an initial error does not affect the order of numerical 
accuracy. 

By simple manipulations we obtain the following equivalent formulation [2]: 

1 + ^ /+i _ sA h ^ k+l / 2 = (l + ^<f> k + 2(3ip k , (2.8) 

= ^ + - (0 fc+1 - 4> k ) , (2.9) 
s 

which shows that the equations may be decoupled. In fact, we can obtain <fi k+1 first by 
solving (|2.8p and then update tp k+1 using (|2.9p . Clearly the solvability of the scheme 
rests on the solvability of Eq. 



2.2. Mass Conservation, Unique Solvability and Unconditional Energy 
Stability. Mass conservation, unconditional unique solvability, and unconditional 
psuedo-energy stability were established in [2]. We recall these facts here, though 
the reader is directed to the reference for details for the details. There are two 
modifications below from what is in [2]. First, our non-dimensional scaling of (jl.ip is 
slightly different, and, second, we use different initializations for our multistep, convex 
splitting scheme. 

Theorem 2.1. The second order MPFC scheme \2. 8\) - \2. 9\) is uniquely solvable 
for any time step-size s > and, moreover, solutions are mass-conservative, i.e, 
U k \\l) = U°\\l), for all k = 1,2,.... 
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Before we state the next result, which is proved in [2], we introduce a third fully 
discrete energy: for each time step k > 1, set 



ik ik-l 



(2.10) 



Theorem 2.2. The second order MPFC scheme h2.8\) - (2. 9\) (or equivalently 
\2. 6]J ) is unconditionally energy stable. In particular, suppose that (f> k , ip k , 
_1 G Cmxn are periodic, and that <fi k+1 , /i k+1 ^ 2 , ip k+1 G Cmxn is a periodic solution 
triple to P3J) - IEM- Th en, for any k>0, 



where 



6 fe+1 ,0 fc ,V fc+1 



fe+1/2 



1 + Yll V fc(W)ll2= J7 (^r '^). 



(2.11) 



D 



2 ±k 



2<j) k 



(2.12) 



This next result follows by summing Eq. (|2.1ip from k = to k = £ — 1. 
Corollary 2.3. With the same assumptions as in Thm. [2~B we have 



fe=0 



1 ±k\ 



k=0 



F(4>°). (2.13) 



Using Lems. [AT41 and IA771 we find 

Lemma 2.4. Suppose that <ft G Cmxn is periodic. Then the following estimates 
hold: 

F(<t>)>C 5 \\4,\\l 2 -^, (2.14) 

m>CkU\t-^ t ^:=g, (2.15) 

m>C r \\V h <i>\\l-^, Cr:=g, (2.16) 

where C5 > and on/y depends upon a. 

Using the last two results and the simple estimate 

F(4> k ) <F{<t> k A h -\Tp h ), (2-17) 

for any > 1, we obtain 

Theorem 2.5. Let <£> be a sufficiently regular, periodic solution to (|1.1| on 
fl T = (0, L x ) x (0 x L y ) x (0, T), with 9 t $( • , • ,0) = and ^ = ^r 1 := $(p i)Pj -, 0), 
■0° = 0. Suppose E is the continuous energy \1.2\) and F is the discrete energy \2.1\) . 
Let <f) k j £ Cmxn be the k th periodic solution of \2. 8\) and i2.9]) for 1 < k < I. Set 

M :=E($(-,-,0))+C s L x L y , (2.18) 
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where Cg > is a constant that does not depend on either s or h. Then we have the 
following estimates: 



e-i 



fc=0 



s 

, + ~2 



k=0 



J2\\V h (D*<t> k )\\l<M , 

lljLfcll / 

max \\4> < a/tt- 



max 

0<fc<£ 



max V;,( 

o<k<e 11 



< 



< 



5 

Mo 




C 6 
Mb 




Cio, 



(2.19) 

(2.20) 
(2.21) 
(2.22) 



Theorem 2.6. Suppose that &(x, y, t) is a periodic solution of the MPFC equation 
jl-4]) , with the regularity assumed in Thm. \3.3\ below, such that <9 f $(x, y, 0) = 0. Then 
we have the following estimates: 



L°°(0,T;L°°(n)) 



ll*lli-(o,r iW ri.«( n)) < W Cie (e ($ (*, y, 0)) + ^) 



C13 
Cl5 
Cl7 



(2.23) 
(2.24) 
(2.25) 



/or any T > 0, where C12, C14, Ci6 > are constants that are independent of T . 

3. Error Estimate for the Second Order Scheme. We now prove an error 
estimate for the second order scheme (|2.8I) - (|2.9p for the MPFC equation. The 
following estimate, proved in |22j . shows control of the backward diffusion term. 

Lemma 3.1. Suppose that (f> 6 Cjaxn is periodic and that Ah<fi G Cmxn is also 
periodic. Then 



\±h<t>\\l < 7^ 



\ 2 2 + ^\\V h (A h < 



| 2 
2 ■ 



(3-1) 



valid for arbitrary e > 0. 

In addition, a control of the error related to the nonlinear term in the second 
order scheme is needed. 

Lemma 3.2. Suppose <£> fe , $ fc+1 , <j) k+1 g Cjnxn are periodic and denote their 
ferences by cf> k := $ fc — cf> k and 4> k+1 := $ fe+1 — tfi k+1 . Then we have 



\A h ( X ($ k+1 ,* k )-x{<l> k+1 ,<t> k ))\ 



<C^Kl- ( 


A h ~4> k+1 


+ 


A h 4> k 


2 )+ k M 


v h 4> k+1 


+ 


v h 4> k 




2 








4 





( 


4> k+1 


+ 


4> k 


00^ 






00 







{Kl + K X K 2 ) • ( 



4> k+1 


+ 


4> k 


J} 




2 







(3.2) 
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with 



Ki = 
K 2 = 

K 3 = 
K 4 = 
K 5 = 



$ k+1 \\ + 

1 1 oo 

A?* fe+1 || 

'i Hoc 

A^ fe+1 || 2 



| <£ fc || _j_ | 



1 1 OO 

|AS$ fc || 

I '* Moo 
A X lk 1 1 I 



+■ Wh® k+1 \ 
\K4> k+1 \\ 2 - 



+ \\K* k \ 

l^lla. 



v fc * 



fe+1 1 



fe+1 1 



|v h $ fc | 
llv h $ A 



|Vfc. 



Js+ll 



(3.3) 



and Ci8 is a positive constant that is independent of h. 

Proof. First, careful expansions yeild the following nonlinear error decomposi- 
tions: 

($ fc+1 ) 3 - (0 fe+1 ) 3 = (($ fc+1 ) 2 + $ fe +V fe+1 + (0 fe+1 ) 2 ) 4> k+1 , (3.4) 

($fe+i) 2 $fe _ (^+i) 2 ^ = + 0*+i) + (^+i) 2 ^ , (3.5) 

($*) 2 _ (/) 2 = ($* + + (^) 2 , (3.6) 

($ fc ) 3 - (0 fe ) 3 = (V) 2 + $ fe fe + (0 fe ) 2 ) 0* . (3.7) 

Meanwhile, a detailed calculation yields the following finite difference expansion: 

(D x gi + i/ 2 jD x hi+i/ 2t j + D x gi_i/ 2 jD x hi_i/ 2 j^ 

+9i,j {D x fi+i/2,jDxhiJ r i/ 2 ,j + Dxfi-i/2,jD x hi-i/ 2 j) 
+hi+i,jDxfi+i/2,jD x g i+ i/ 2 j + hi-ijD x fi_i/ 2 jD x gi_i/ 2 j . (3.8) 

An analogous formula for A^ (fgh) i j holds by symmetry. First, we bound all of the 
terms in the expansion of (($) 3 — (</>) 3 )- For brevity, we only show how this is 
done for one term, namely, A£ (^((f> k ) 2 4> k ^j ■ The expansion is given by 



A 



N^+2nQ+2nQ+N^, 



with 



(3) 



(2) 



N. 



1,3 

(4) 
i,3 



lj (D x 4> k +1/2tJ Dj k +1/2 . 3 4 D 



k i+i,j 



D 



V-1,3 



x<Pi-l/2,j D x<f>i-l/2,j) 
k 



D x c, 



(3.9) 

(3.10) 
(3.H) 
(3.12) 



Discrete Holder's inequalities can be applied to bound all of the above terms as follows: 



N {i) 

N (2) 
N (3) 



< 



< 



< 2 



< 



<2 V h <f> k 



(3.13) 
(3.14) 
(3.15) 
(3.16) 
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with repeated application of Lem. IA.51 The nonlinear error term (($) 3 — (0) 3 ) 
can be analyzed in exactly the same way. Combining the estimates using the triangle 
inequality gives the result (|3.2I) and the lemma is proven. □ 

We now establish an error estimate for the fully discrete second order convex 
splitting scheme for the MPFC equation. We do this in three steps. First, we derive a 
local truncation error for a finite Fourier projection of the exact solution to the MPFC 
equation (jl.ll) . Second, we derive an estimate of the difference between our numerical 
solution to the scheme (|2.4j) - (|2.6j) and this finite Fourier projection. Third, we use 
the triangle inequality to derive our global error estimate. 

In the rest of the paper, for notational simplicity only, we will assume L x = L y = 
L, and hence SI = (0, L) 2 . As a consequence we have m — n — N , where we may as- 
sume N is even. The more general rectangular case can be handled straightforwardly. 
Now, suppose that $ has the following Fourier series representation on SI: 

oo 

*(x,y,t) = ®kAt)^ {kx+lv) , (3.17) 

k,l= — oo 

with 

=p / ^(x,y,t)e-^ kx+l ^dxdy . (3.18) 

The (finite Fourier) projection of $ onto the space B N I 2 , consisting of all trigonometric 
polynomials in x and y of degree up to N/2, is defined as 

N/2 

® N {x,y,t):=V N ®[x,y,t):= ]T $ M (i)e^ fe *+^ . (3.19) 

k,l=-N/2+l 

Define 

s 2 

V N (x, y, t) := d t $ N (x, y, t) - —df$ N (x, y, t) . (3.20) 

For any function G = G(x,y,t), given s > and k > 0, we define G k (x,y) :— 
G(x,y,s- k). 

Theorem 3.3. Suppose the unique periodic solution for the MPFC equation {1.4]) 
is given by 

$ £ H 4 (0, T; L 2 (SI)) n L°° (0, T; H 8 (SI)) 

HW 2 ' 00 (0,T;H 2 (Si)) C\H 2 (0,T;H 6 (Si)) , (3.21) 

for T < cxd. Set W := d t <£>. Then 

(3.22) 
(3.23) 



/v 



2 V 



fe+i 

JV 



4 A ' 



fc+i 



AT 



1 



fc+1 



N 



A h (3$ 
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where T k and satisfy 

\\ T i\\L 2 a (Q,T;L 2 h (fl)) 



T/s 
\ k=0 



„k+l I 



< M (s 2 + h 2 ) 



(3.24) 



for i = 1,2, with 

M < C(\\3>\\ H 4( 0T . L 2(Q^ + \\®\\w 2 -°°{0,T;H 2 (n)) + ll^llff 3 (0,T;H 6 (n)) 



1*1 



L°°(0,T;ff 4 (O)) 



(l + 11*11 



H 2 (0,T;H 2 (Q,)) 



L°°(0,T;ff 8 (Q)) 



(3.25) 



The details of the proof are technical and are contained in the appendix. 
Remark 3.4. The constructed solution \3. 1 9\) comes from the Fourier projection 
of the exact solution $. The reason for the choice of $n instead of <fr is the fact that 
<S> N G B N ' 2 , which in turn gives a local truncation error estimate without involving 
with an aliasing error, as can be seen in the appendix. Meanwhile, an 0(s 2 ) correction 
term is added in the construction LS.20\) for so that a higher order consistency is 
obtained in A3. 23\) . Such a correction term is based on an asymptotic expansion of the 
numerical scheme and the resulting higher order consistency is crucial in the stability 
and convergence analysis. Finally, a numerical convergence of the numerical solution 
to *&n) is equivalent to its convergence to the exact solution (<£>, *), since $ is a 
spectrally accurate approximation, and 'J n is an 0(s 2 ) approximation to "J*. 

Theorem 3.5. Suppose $, $jv ; *B and ^ n are as in the last theorem. Define 
<t>li ■= ® k N (h-i,h- j)-<l>i,j and $v := V k N (h-i,h- where ^ G C„xn 

are the k th periodic solutions of \2.J$ - (2. 6\) . or equivalently, \2. 8\) - &2.9\) . with 
<fq j and ipf j = 0. Then 





"*,3 



<C(s 2 + h 2 ) , 



(3.26) 



provided s is sufficiently small, for some C > that is independent of h and s. 
Proof. Subtracting ([23) - from ([3T22]) , (j3~^3l) yields 



1> 



k+1 



l\) k 



lk+1 



= A h \ [ x{® N ,* N 



hk+1 



k-1 



A^+1 _ , 



(3.27) 
(3.28) 



where </> fc+1 / 2 := \ (^> k+1 + ^ and ^ k+1 / 2 := \ ^jj k+1 + . Taking the inner 
product with the error difference function h 2 (j> k+1 — (f> k ^j gives 

h 2 ^-^ T k )+h 2 (0*+i_0*||A fc ( X (**+ 1 J ^)- X (^ 1 ,0*))) 



2* U k +i p 

s \ 

-ah 2 (> +1 - 4> h 



Ik+i 



— U^ 1 -^ 

s \ 

- Uk+l _ ~^k 



lk+1 



A 3^+l/2 



-h* cp* 



lk+1 



k-l 



(3.29) 
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The first term on the right-hand-side of (|3.29l) can be rewritten and estimated as 
follows. With the help of (|3.28l) and an application of Cauchy's inequality we have 



s 

_ 1 
~~ 2 
1 

> - 
~ 2 



(> +1 - ft 



4> 



k+l 



■0 fe+1 + ip k 



J + sh 2 (7 



2 S H T 2 II a - s 



2 



k+l 



V 7 



The second term on the right-hand-side of (|3.29[) is obviously non-negative: 



s 



(> +1 - ft 





ft +1 - ft 


ft +1 - ft) = s 




s 



> . 



(3.30) 



(3.31) 



The first term on the left-hand-side of (|3.29l) can be controlled using Cauchy's in- 
equality: 



h 2 u 



lk+1 



D = *h 2 



fc+1 



1j) k 



< 



k+l 



'1 II 2 



'2 || 2 



.(3.32) 



The analysis of the convex diffusion terms can be carried out with the help of the 
discrete Green's identities (|A.8I) and (| A.9|) : 



lk+1 



A h ft +1 ' 2 



k+l 



and 



(> +1 - ft 



A 3jk+l/2 



V h (A h ft +1 ) 



V h (A h ft 



The concave diffusion term can be handled with the identity 

2 1 



h 2 



lk+1 



Ife— 1 



A h 6 k+1 



A h ft 



2 2 
2 1 



lk+1 



(3.33) 



(3.34) 



A h (ft -ft- 1 ) 



A 



h (ft +1 - 2ft + ft- 1 )^ . (3.35) 



For the nonlinear term, we start with an application of Cauchy's inequality: 



Tfc+l _ Tfe 



= 8* 



< 



k+l 



k+l 



A h { x ^ k N +\^ k N )- x (ft+\ft)} 



\A h { x ^ k N + \^ k N )- x (ft+\ft)}\ 



(3.36) 
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Lem. 13.21 can be used to bound the last term appearing above. In more details, the 
following uniform (in time) estimates are recalled from our previous lemmas: 



N\ 



oo 



<||*y L oo<Cl6 



< c v$ 



n\\l° 



' •' !: - r ' " • " <CC 9 , ||Ajy|| 2 <C 
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C < C , (3.37) 
< C 9 , (3.38) 
(3.39) 



and the following estimates are valid on the finite time interval [0, T] 



< ||w|| ioo + c< c , 



A?$' + A?*' <\\d xx $ 

11 Moo I " M oo — M J " i 



c <c 



(3.40) 
(3.41) 



for I = k, k + 1, where C denotes a generic positive constant that is independent of h. 
Applying Lem. IA.7| Lem. IA.4| and substituting estimates (|3.37[) - p.41|) yield 



A h { X ($ k N + \$ k N )-x(<f> k+ \<f> k )} 



< C 



19 





+ 




+ 




+ 




J 




2 




2 




2 







(3.42) 



where C19 > is independent of h and s, but is dependent upon T and also the exact 
solution $. Going back to (|3.36[) and using the last estimate and Lem. 13.11 (with 
e = 1) we obtain an estimate for the nonlinear term: 



(> +1 



A^xK+^-x^ 1 ,/)} 



s 

< - 
~ 2 



fe+i 



■2s 2 IN 



-2s C 



19 



( 




2 

+ 


fc 


2 

+ 


A,^' +1 


2 

+ 


A,# 








2 




2 




2 







< 



fe+i 



■2s 2 IN 



V h A^ fe+1 



(3.43) 



Define a modified energy for the error function via 



o 



2 2 

fe Ifc-i 



1 



2 2 
2 



V, 



A combination of ([3~29)) . ([3"3Tj) - (|3~35|) and (pT43j) results in 



(3.44) 



(4> k+1 )-F 1 ^ k )<sC 2 



+ 



fc+1 



V h ( A^ fc+1 



(A^ fc ) 



'2 lb 



(3.45) 
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where C 2 q > is independent of h and s. Summing over k and using the fact that 



Fi 



< Ch 4 yields 



Fi (>) <2 S C 20 ^ 



k=l 



+ 



h (a,^) 



'•- ( ' ( II T iIIl2( T;L 2 (o) ) + ll T 2 11^2(0 T;L 2 (0) ) 



< 



+ 



k=l 

+CM 2 T{s 2 + h 2 ) 2 . 



(3.46) 



To carry out further analysis, we introduce the positive part -Fi: 



F- 



V h cf> k 



2 1 

+ 2 



V h A 



A, 



(3.47) 



so that (|3.46|) becomes 



< 



k=l 



7 h (AhP) 



■„2\2 



+ \\A h </> e \\ + CM 2 T{s 2 + h 2 ) 



To estimate the additional term 
The following identity is observed: 

Ik Ik-i 



we need a bound of 



(3.48) 
in terms of ip k . 



k=l 



= d>" + s 



E 

fe=i 



k-l 



'2 . 



(3.49) 



with error equation (|3.28[) used in the last step. In turn, an application of Cauchy 
inequality shows that 



< 2 



-4sT^] ||$* 



fc=i 



4 S J T||r 2 | 



L2(0,T;Lj(n)) 



< 4sT^ |U fe + C(/i 4 + s 2 ( s 4 + /i 4 )T), 



k=l 



(3.50) 



in which the fact that 0° < C/1 2 (which comes from the construction (13.191) of 
II II2 

the approximate solution and the initial numerical data j = along with the 

truncation error analysis (|3.24j) . Therefore, using Lem. 137X1 shows that 

,||2 2e 



2 1 
2 " 3? 



(a^) 



< 



2s T 
"3?" 



fc=i 



2 2e 
2 3 



+ C(s 4 + h 4 ) , (3.51) 
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for any e > 0, with a trivial requirement that s 2 T < 1. Taking e = |, the substitution 
of the estimate (I3.51|) into (I3.48|) shows that 



1 2 1 



+ 



k=l 

+CM 2 T(s 2 + h 2 ) 2 , (3.52) 
where C21 > is independent of h and s. Introducing the more refined energy 



'h (A^") 



2 2 



1 

+ 2 



2 1 

2 + 4 

2 



V,,. (A^ fc ) 



= F 2 U K -- 



Vh(A^*)| , (3.53) 



we obtain, with the aid of the estimate (|3.50[) . 



t k 



F 3 < sC 22 ^f 3 (4> k ) + s 2 T C 21 Y,Y, \ +CM 2 T(s 2 + h 2 ) 2 , (3.54) 

fc=i k=ii'=i 2 

where C22 > is independent of h and s. Meanwhile, motivated by the estimate 

£ k £ k £ 



k=l i' = \ 

which follows from 



k=l l' = \ 



k=l 



' = 1 



E ¥' , * E ¥' 



'=1 



, V k < I 



we arrive at 



F: 



(?) 



k=l 



CM 2 T(s 2 + h 2 ) 2 



(3.55) 



(3.56) 



(3.57) 



where C23 > is independent of h and s. Applying a discrete Grownwall inequality 
gives 



F 



s(^) <C 24 (.s 2 + ^ 2 ) : 



(3.58) 



which holds provided s is sufficiently small. Note that C24 is a positive constant that 
is dependent upon T (exponentially) and but is independent of h and s. □ 

Corollary 3.6. Define := <fr fc (pi,pj,ks) - </>?„.. T/ien 



+ 



V,, 



(A^ fc )|| 2 <C( s 2 + / i 2 ) , 



(3.59) 



provided s is sufficiently small, for some C > that is independent of h and s. 

Proof. Estimate (I3.58|) gives the discrete H 3 estimate for <\>. For the projection 
$n, we have the following approximation estimate: 



\^N~n L ^(0,T,H" 



< Ch m ||$| 



L°°(0,T;H 



(3.60) 
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for m, r > 0. See for example the references [3] [3] [IT]. Combining estimate (|3 . 58[) 
with the approximation result p. 601) . we can obtain the estimate (|3.59p . □ 

Remark 3.7. By virtue of h3. 59}) and Lems. \^4\ and [3~l[ along with the estimate 
Ii3.50\) , we immediately get an error estimate of the form 



<C(h 2 + s 2 ) 



(3.61) 



4. Conclusions. In this paper, we have established the convergence analysis of 
an unconditionally energy stable second order accurate finite difference scheme for 
the sixth-order Modified Phase Field Crystal (MPFC) equation. The parabolic Phase 
Field Crystal (PFC) equation, which is a mass conserving gradient flow, is obtained 
as a special case of the MPFC equation. The numerical scheme is based on a second 
order convex splitting of a discrete psuedo-energy and is semi-implicit. 

Appendix A. Tools for Cell-Centered Finite Differences. 

In this first appendix, we define the summation-by-parts formulae, discrete norms, 
and estimates in two space dimensions that are used to define and analyze our finite 
difference scheme. With some exceptions, the theory will extend straightforwardly 
to three-dimensions. Here we use the same notation and results for 2D cell-centered 
functions as from Sec. 2]. The reader is directed there for all of the missing 
details. 

For simplicity, we assume that f2 = (0,L X ) x (0,L y ). The framework that we 
describe has a straightforward extension to three space dimensions. Here we use 
the notation and results for cell-centered functions from (2TJ [22]; see also [10[ [19] . 
The reader is directed to those references for more complete details. We begin 
with definitions of grid functions and difference operators needed for our discretiza- 
tion of two-dimensional space. Let Q = (0,L X ) x (0,L y ), with L x — m ■ h and 
L y — n ■ h, where m and n are positive integers and ft, > is the spatial step size. 
Define p r := (r — l /2) ■ h, where r takes on integer and half-integer values. For 
any positive integer £, define Eg = \p r I r = |, . . . ,1 + 5}, Ct = {p r \ r = 1, . . . ,£}, 
Cj = {p r ■ h I r = 0, ...,£ + 1}. Define the function spaces 

Cmxn = {4> ■ C m x C n -)• R} , C- lXlT = {</> : C m x -> K} , (A.l) 

Cmxn = {<P : Cjjt x C n -> R} , C myM — {<$ : C m xCj-)R}, (A.2) 

£™ xn = {u:E m x C n -> R} , £™ xn = {v : C m x E n -> R} , (A.3) 

Cxi = {u : E m x C w -> K} , £^ xn = {v:C m X E n -> R} . (A.4) 

We use the notation <fii t j := <fi(pi,Pj) for cell-centered functions, those in the spaces 
Cmxn, Cmxn, C m xn, or Cmxn- In component form east-west edge-centered functions, 
those in the spaces £™xn or ^mxni are identified via Ui + y 2 j := u(p i+ i/ 2 ,pj). In 
component form north-south edge-centered functions, those in the spaces £mxm or 
£^ xn , are identified via u i+ i/ 2 j := u(pi + i/ 2 ,Pj). The functions of V m xn are called 
vertex- centered functions. 

We need the weighted 2D grid inner- products ( • || • ), [ • || • ] ew , [ • || • ] ns that are 
defined in [21] [22]. In addition to these, we also need the following one-dimensional 
inner-products: 

m n 

(/*J + l/2 |.9*.j + l/2) = fj,j+y29i,j+ 1 /21 (/i+l/2,*|<7i+y2,*) = fi + l/2,j9i + l/2,jl 

i=l j=l 

(A.5) 
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where the first is defined for f,g£ £° s xn , and the second for /, g £ £™xn- 

The reader is referred to [3TJ [35] for the precise definitions of the edge-to-center 
difference operators d x : £™ xn — > C mxn and d y : £JJf xn — >■ C mxn ; the x— dimension 
center-to-edge average and difference operators, respectively, A x , D x : Cmxn — > ^mxrv 
the y— dimension center-to-edge average and difference operators, respectively, A y , D y 
Cmxn — > £mxni an d the standard 2D discrete Laplacian, A^ : Cmxn — > C m xn- 
The summation by parts formula we need from |22] are the following: 
PROPOSITION A.l. (Summation-By-Parts:) If <fi £ Cmxn UCm X5J and f £ £%? xn 
then 

h 2 [A^H/U = ~h 2 (<P¥ x f) 

-h (A x (j)x/ 2 ^\fi/ 2 ^) + h (A x (f>m+i/2,i,\f m +i/2,*) , (A. 6) 

and if 4> £ C m xn U Cmxn and f £ £ r " s x „ then 

h 2 [D y m ns = -h 2 mdyf) 

-h (Ay^i/^f^i/A + h (A y <j)*m +1/2 \f+^ n+1/2 ) . (A.7) 
Proposition A. 2. (Discrete Green's Identities:) Let (f>, tp £ Cmxn- Then 

= -h 2 (<j>\\A h tp) - h (A x (/)i/ 2t± \D x ipi /2jic ) + h (A x (f> m +i/2,*\Dxif>m+i/2,±) 

-h {Ay^x/^Dy^i/A + h (A v ^ )n+ y 2 |D y V*,n+Va) ' ( A - 8 ) 

and 

h 2 (0||A^) =h 2 (A^||^) 

+h (A x (j)m+i/2,*\D x 1pm+-i-/2,*) —h (D x (f)m+i/2,*\A x 1pm+i/2,*) 

-h (A x c/)i/ 2 ^\D x i(!i /2 ^) +h (Dx<l>ya,*\Axi>i/i,*) 

+h (A y (/)^ yn+ if 2 \D v ip^ 7n+ i/ 2 ) -h (Dy^^+i/^Ay^^ 

-h (Ay(f>*i/ 2 \D y ^i/ 2 ) + h (Dytp^i/^Ay-ip+i/z) . (A. 9) 

In this paper we are interested in periodic grid functions. Specifically, we shall 
say the cell-centered function (f> £ Cmxn is periodic if and only if 

(f>m+i,j = <t>i,j, <t>o,j = <t>m,j> 3 = 1) • • • (A.10) 
4>i,n+l = 0i,l' 0j,O = <f>i,n, i = 0,...,m + l. (A. 11) 

For such functions, the center-to-edge averages and differences are periodic. For ex- 
ample, if £ Cm X n is periodic, then A x (f> m+ i/ 2 j = A x (f>i/ 2 j and also D x (f> m+ x/ 2 j = 
D x <f>y 2 j, for all j = 0, 1, ... , n+1. We note that the results for periodic functions that 
are to follow will also hold, in a possibly slightly modified form, when the boundary 
conditions are taken to be homogeneous Neumann. 

We define the following norms for cell-centered functions. If £ C m xn, then 
||0|| 2 '■= \/h 2 (0H0), and we define |]V/i0|| 2 , where 6 Cmxn, to mean 

l|V h 0|| 2 := ^ 2 [A^HA^U + Zi 2 [A,0||A,0] ns ■ (A.12) 
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We will use the following discrete Sobolev-type norms for grid functions <f> G Cy, 

\\<l>\\o,2 : = Hh and 



H\\i,2--=\/U\\ 2 2 + Wf 2 , U\\ 2 . 2 --= y/nWl + WVhtWl + lMWl ■ (A.13) 



In addition, we introduce the following discrete L and L°° norms: for any <fi G C n 
define 



{h 2 (m) 



1/4 



and 



(A.14) 



And for G Cmxn define 

||V^|| 4 := (h 2 [(D x cj>) 4 



1/4 



(A.15) 



Some discrete Sobolev-type inequalities for two-dimensional grid functions are 
needed in the analysis in later sections. The following results are recalled; the detailed 
proofs can be found in [22l [19] . 

Lemma A. 3. Suppose that <f> G Cmxn- Then, 



W 4 <CiWi, 21 Ci:= 2max 



max <J — . L , } . max <{ J^i L y 



1/4 



(A.16) 



Lemma A. 4. Suppose that <p G Cmxn is periodic. Then, for any i G {1, 2, . . . , m} 
and any j G {1, 2, . . . , n}, 



|^| 2 <C 2 ||< 2 , C 2 :=4max{-i-,^-,^,^ 

Hence U\\l<C 2 U\\l 2 . 

Lemma A. 5. Suppose that <f> G Cmxn * s periodic. Define 



(A.17) 



where 



S ■= h " E E «' D v (^Wy+y* 

i'=0 j'=0 



1 if k G {1,2,...,*- 1} 

1/2 ;/ fcG{o,n 



(A.18) 



(A.19) 



Then S — h 2 (Af^W A?^), where := d x D x and A v h := d y D y are the 3-point discrete 
lapacian operators in the x- and y-directions, respectively [22]. And, since S > 0, we 
have 

h 2 (Al<f>\\A%cf>) < h 2 (A^HA^) and h 2 (A^||A^) < h 2 (A h ^\\A h cP) . (A.20) 



Consider the space 



H:= {<peC mxn \(<j>\\l) = 0}, 



(A.21) 



SECOND ORDER SCHEME FOR THE MPFC EQUATION 



17 



and equip this space with the bilinear form 

for any (f>i, <p2 G H, where ipi G Cmxn is the unique solution to 
- A h ipi =(f>i, V* periodic, (ipi\\l) = 0. 



(A.22) 



(A.23) 



Lemma A. 6. The bilinear form (0x i * s 071 inner product on the space H. 

Moreover, 



{M<hU = - (^illA^ 1 (&)) = - (A^ 1 (<t>i)\\fo) 



Thus 



(A.24) 
(A.25) 



defines a norm on H . 

Lemma A. 7. Suppose that <fi G Cmxn "is periodic and set 

||0 - 0|L < c 3 IIV^^IL , 



— i — (0111). Then 

m-n Vr II / 



(A.26) 



where C3 > is a constant that only depends upon L x and L y . Furthermore, 

||0-^|| 4 <Ci||V^|| 2 , ||V^|| 4 <Gi||A/,^|| a . (A.27) 

where C 4 := Ci^/CfTl. 

Appendix B. Consistency Analysis of the Second Order Numerical 
Scheme. 

In this appendix we give a detailed derivation the local truncation error estimate 
(I3.25[) . We establish the results for vertex-centered grid functions, rather than cell- 
centered functions, as the indexing becomes simpler. 

B.l. Proof of Estimate (|3.25[) . The following three results will be used to 
establish (ET25)) . 

PROPOSITION B.l. For f G H 3 (0,T), we have 



t rk 



T f\\mo,T) ^ CsTn ll/llff^+w) . mth T f 



f k+i - f 



-f{t k+1 i 2 ) 



(B.l) 



for < m < 2, where C only depends on T, \\ ■ ||x,2(qx) * s a discrete L 2 norm (in 

time) given by ||s|| L 2 (0jT) = (9 k f ■ 

PROPOSITION B.2. For f G H 2 (0,T), we have 



Dt /2 f 



L- 



[T/s]-l 2 

« E ( D t/2f k+1/2 ) <c\\f\\ HH0 . T) 



I D 2 f II 

I t J \\l%(o,t) 



(0-T) ' \ 

1/2 _ 4(/ fc + 1 -2/(.,t fc +V 2 ) + / fc ) 



(B.2) 



with D 2 t/2 f k 



\ 



[T/s]-l 

s (D 2 f k f<C\\f\\ 

k=0 



H 2 (0,T) 



(B.3) 



with D 2 f = 



2 fk / fc + 1 -2/ fe + / fe - 1 



18 



A. BASKARAN, J.S. LOWENGRUB, C. WANG, AND S.M. WISE 



where C only depends on T . 

The proof of Proposition IB.1IIB.2I is based on the integral form of the Taylor 
expansion in time. The details are skipped for the sake of brevity. In the spatial 
discretization, the following proposition gives a corresponding 0{h 2 ) truncation error 
bound. 

Proposition B.3. If f e B N I 2 has a regularity f e H^ er (fl), we have 

\\A k f-A h h f\\ Ll{n) <Ch 2 \\f\\ H2+2km , for k = 1,2,3, (B.4) 



where C only depends on Lq and WgW^^ = >J h? X^ij=o 9i j • 

The key point of this proposition is that the projection approximation solution 
$at S B N / 2 so that an aliasing error is avoided in its centered difference approxima- 
tion. This consistency analysis can be carried out by a detailed Fourier expansion of 
both A^$jv an d A$jv at the discrete level, and the comparison of the corresponding 
discrete Fourier coefficients leads to the above estimates, following a similar method- 
ology as in [T5]. The details are skipped for brevity. 

Observe that the O (s 2 ) correction in the definition of (|3.20|) is added so that a 
higher order consistency between at £ fe+1 / 2 anc j (cf)^ 1 — $^r)/s can be derived. 
Looking at the time derivative of the projection operator, we observe that 

Qk Qk 8 k $(x t) 

In other words, 9*$jv is the truncation of 9 t fc $ for any k > 0, since projection and 
differentiation commute. This in turn implies an accurate approximation of the cor- 
responding temporal derivative, at any fixed time: 

\\d k (3> N - < Ch m \\d^\\ Hm+r , (B.6) 

for m, r > 0, and < k < 2. 

Since the exact solution of the MPFC equation has the regularity (|3 . 2 1 1) . the 
approximation estimates (|3.60I) - (|B.6I) imply the same regularity for the projection 
solution $jv, vffjv, by taking m = 0: 

\\$N\\ H r <C ll^^-^II^^Cll^H^, (B.7) 

at any fixed time. 

We define the following quantities: 

fc+l/2 _ * k r + 1 -$>% 



Fx 

^fc+l/2 
b 2 
? k+l/2 



k + 1 / 2 _ A3(b Hl/2 

'h N 
fc+l/2 _ gg^+frjy 



^fc+l/2 
b Xe 


= {9?$n) ( 


■ ,t k+1 / 2 ) , 


pfe+1/2 
b 2e 


= (d t $ N ) ( • 


>t *+l/2) ( 


vk+1/2 
*3e 


= A ((<&at) 3 ) 


( • ,i fe +V2) 


jpk+1/2 


= A(*jv)( ■ 


^fc+l/2) f 


pfc+1/2 


= (A 2 $ w )( 


. jt fc+l/ 2)) 


^fc+l/2 
P 6e 


= (A 3 <S> N )( 





F 3 fe+1/2 = A, (($ 3 N ) k+1 ? 2 ) , 

F 4 fc+1 / 2 = A, (^ +1 / 2 ) , F^ 2 = A(**) ( ■ ,i*+V») , (B.S 

F k+ ^ = A 2 (f 4* - i^ 1 ) 



^7—2 

Moreover, the corresponding values for the exact solution are denoted by 

F k ;tn' 2 = dM ■ ,^' +1/2 ), f&i* = om ■ 

= (A$ 3 ) ( • ^+1/2 = (A$) ( . )t fc+l/2 }) 

F^ 1 ' 2 = (A 2 $) ( • ,«*■*/»), ^+1/2 = (A 3 $) ( . ^+1/2). (B _ 9) 
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Note that all these quantities are defined on the numerical grid (in space) point-wise. 

First we look at the first order time derivative term, F 2 , F 2e and F 2en . A direct 
application of Proposition IB . 1 1 indicates that (by taking m = 2): 



F 2 - F 2e \\ mQtT) < Cs 2 \\^n\\ h3(0!T) < Cs 2 \\<S>\\hHo.t) - 



3 (0,T) 



(B.10) 



for each fixed grid point (i, j), in which the second part of (|B.7I) was used in the second 
step. Meanwhile, the approximation estimate (|B.6|) yields (with k = 1, m = 2): 



F, 



fc+1/2 



2c 



F. 



fc+1/2 



2en 



< Ch 2 \\d t $\\ H2 



(B.ll) 



Therefore, a careful calculation shows that a combination of the above two estimates 
results in 



\Fo - F, 



2en\\Ll(p,T;Ll(Sl)) 



<c( s 2 +h 2 ) (\m 



ff 3 (0,T;L 2 ) 



1*1 



(B.12) 



Similar analysis can be applied to the second order time derivative terms. The 
construction Q3.20p for the approximate solution gives 



F, 



fe+l/2 



9 t $t +1 - d t <5> k 



N 



s 2 d?$ N 



fc+1 



N 



12 



F, 



fc+l/2 



pfe+1/2 



12 



(B.13) 



in which Fu and Fi 2 are finite difference (in time) approximation to df$>N, df^N, 
respectively. In more detail, if we denote F 1 ^^ 2 — 8 2 ^n( • ,t k+1 / 2 ), f\ 2 ^ 2 = 
df$N{ ■ ,t k+1 l 2 ), the following estimates can be derived by using Proposition IB. II 
(with m = 2 and m = 0): 

11*11 ~ *iielL| (0 ,T) < Cs 2 \\* N \\ H * %T) < Cs 2 ||$|| ff 4 (0 , T) , (B.14) 

11*12 - *12e|li3( ,T) ^ C\\$n\\h*(0,T) - C H $ Hff 4 (0,T) ' (B.15) 

for each fixed grid point Again, the approximation estimate (|B.6I) gives (with 

k = 2, m = 2): 



pfc+l/2 _ pfe+1/2 



<Ch 2 \\d 2 $\\ H2 . 



(B.16) 



A combination of (|B.13|) - (IB.16|) leads to 

||*i - F len \\ L2 ^ T . Ll{n)) < C(s 2 + h 2 ) (||$|| ff 4 (0iT;i2) + ||$|| M ,2,oc (0)T .^ ) )B.17) 

For the convex diffusion term F4, F± e and Fi en , we start from an application of 

fc+i_ 



Prop. ES (recall that <S> k + 1/2 = 2 +$JV 



F 



fe+l/2 



(^ +V2 ) 



< Ch 2 



fe+l/2 
N 



< C/l 2 ||^iv|lioo( 0;r . H 4) 



Meanwhile, a comparison between ^j^ - ' and $n{ ■ ,t /c+1 ^ 2 ) shows that 

^ +1/2 -^(-,* fc+1 / 2 ) = ^ 2 A 2 /2^ +1/2 - 



(B.18) 



(B.19) 
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On the other hand, an application of Prop. HOI gives 



N 



<C\\A$ N \\ H2 



(0,T) 



(B.20) 



at each fixed grid As a result of (|B.18[) - (IB.20|1 . we get 

11*4 - F ie\\ L 2(o,T;Ll(n)) ^ + h *) ( II *^ \\ L oo ( ,T;H*) + II II ff 2 (0,T;ff 2 ) 



< + fr 8 ) M|$|| i « ( o, r . H 4) + ||*II^( ,T; 



(B.21) 



The approximation estimate of to i^en is straightforward, from p.60[) (with m = 
2): 

||F 4e - *4en|| L 2( 0T;i 2 (n) ) < Ch 2 II^Hiocjo.TjH*) ■ (B.22) 

Consequently, we arrive at 

||*4 -*4en IIl 2 (o,T;L 2 (Q)) ^ C(« 2 + ^) (||*|| L 

11*11^(0,^)) (B-23) 

The other convex diffusion terms i*6j *6e and i^en can be analyzed in the same 
way. The details are skipped for simplicity. 



11*6 " F 6e n\\ Ll {o >T , Ll{a) ) ^ C ( S + h ) [\\*L-(0,T;H») + 11*11^(0,^) J (B.24) 

The analysis for the concave diffusion terms F§ , F$ e and F^ en is similar to that of 
the convex diffusion term; yet more details are involved. An application of Prop. IB. 31 
gives 



pk+l/2 _ ^2 



k-l 
N 



< Ch 2 



, 1 
$ N ~ -* 



k-l 

N 



H e (n) 



< Ch 2 ||$ 



N\\l°°(0,T;H0) 



(B.25) 



Meanwhile, a comparison between |*^ — 1 and <J>at( • ,t k+1 l 2 ) reveals that 



fe-i 



2 -iv 2 -N 

Similarly, applications of Prop. [BT21 imply 



2 t/2^jv 2 * JV ' 



(B.26) 



D^A 2 ^ 
D 2 A 2 $ r 



< C A^$ 



< C||A 2 $jv 



H 2 (0,T) ' 



L 2 (0,T) 

at each fixed grid Then we obtain 

11*5 ~ *5e|lia(o,T;L=(n)) < C(s 2 + ^) ( II *JV || L « (0,T;ff«) + II ®N II ff 2 (p,T;H*) 
< C(S 2 + h 2 ) (||$|| L oo (0W) + ||*||H»(0,r;H«) 



(B.27) 



(B.28) 



The approximation estimate of J^e to i*5 en can be derived in the same manner: 

||*5e - *5en||i|(o,T;L^(n)) - Cfl ll*IL°°(0,T;.ff<») ■ (B.29) 
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That gives the consistency estimate for the concave diffusion: 

11^5 - F 5en \\ L2s{QT . Ll{n)) < C{s 2 + h 2 ) (\\n L ~{ ,T;HS) + INIi^O.T;^)) ( B - 3 °) 



Next we look at the nonlinear term. A direct application of Prop. IB.3I indicates 



that 



F 



k+l/2 



-A 



((^) fc+1 



/2 



< Ch 2 



3 \fc+l/2 



< Ch 2 (||*^ +1 H 3 ff . (n) + K||*« (n) ) -(B.31) 



in which a product expansion and a Sobolev imbedding are used in the second step. 
Subsequently, we need to compare {^ N ) k+1 / 2 and ($)y-)( • ,t k+1 / 2 ) and derive an 
estimate. A careful calculation reveals that 



($fv)fc+ l/2 _ ($3y)( . ^ + l/ 2) = 1 (( ^ + l )2 + ( ^ )2) . l s 2 {p^ N ) 



k+l/2 



8 



t/2*AF 



k+l/2 



k+l/2\ 



Meanwhile, by the following observation of a nonlinear expansion: 

Mfg h ) = f9 Ah + f hA 9 + <J hA f + 2/V.gV/i + 2.gV/V/i + 2/iV/V.g, 
we obtain 

' ($3^+1/2 _ ($3v)( . )t *+l/2 



(B.32) 



(B.33) 



<CC(C+l)s 



k+l/2 

H 2 



(d 2 /2 <s> 2 n ) 



k+l/2 



H 2 



With C — ||^Jv|li<x>( 0)T . W 2,oo( n )) < C ||$jv||.L°°(0,T;.H- 4 (n)) 

Subsequently, applications of Prop. [B~2l imply 

< C ||*7v||/i-2 (0:T;ff 2) , 



Dt/2®N 



L 2 (0,T;H 2 ) 



D 2 t/2($N) 



L 2 {0,T;H 2 ) - Cf H $ ^llH a (0,T i fl>) 



(B.34) 

(B.35) 
(B.36) 



Note that the second estimate is involved with a nonlinear term <fr 2 N . A detailed 
expansion in its first and second order time derivatives shows that 

8t(& N ) = 2<Z> N d t <Z> N , d 2 {<S> 2 N ) = 2<S> N d 2 <S> N + 2(3 f $jv) 2 , (B.37) 

which in turn leads to 

ll $ ^llff 2 (0,T) - C (ll $ ^Hi°°(0,T) • \\®n\\ H 2 {0 ^ t) + ||$iSr||w-i,4( 0l T)) 

<C\\$ N \\% 2(0>T) , (B.38) 

at each fixed grid point with a 1-D Sobolev imbedding applied at the last step. 

Going back to (|B.36|) gives 



L 2 (0,T;H 2 ) 



< C||$ 



N\\h 2 (0,T;H 2 ) 



(B.39) 
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Therefore, a substitution of (|B.35|) and (|B.39I) into (|B.34I) yields 
A ((<^) fc+1 / 2 -(*%)(. ,t k+1 / 2 )) 



L 2 (0,T;L 2 (fi)) 



<c im 



(II 



•N 



L°°(0,T;ff 4 (n)) 



11*^11^(0^.^2) + l) 



l),s 2 . 



(B.40) 



A' 



For the comparison between F$ e and F^ en , we cannot apply (|3.60[) directly, since $ 
is not in B N / 2 . We observe the difference between and $ 3 is given by 

®n-<S> 3 = (<S>n + + ® 2 ). (B.41) 

As a result, taking a Laplacian operator to the above terms, applying the nonlinear 
expansion (|B.33|) . we arrive at 



pfc+l/2 _ pfe+1/2 
3e 3en 



Z,= (Q) 



<C 11$ 



(ll*A f llL~(0,T;W2 > ~(O)) + II* 



L°°(0,T;W 2 .~(O)) 



• 11$ 



iV 



4>l 



< Ch 2 ||$" 3 



L°°(0,T;ff 4 (n)) 



L~(0,T;ff 2 (f2)) 

(B.42) 



Furthermore, a combination of (|B.31I) . (|B.40[) and (|B.42I) leads to the consistency 
estimate of the nonlinear term 



||-f*3 - ^3 e n|lL2( 0iT;i 2 (f2)) 

< C (s 2 + h 2 ) (||*||io. ( o jT ;jr«(n)) + \\®N\\ 2 L <*, (0 ,T;H*(n)) • ll $ Jvll^(o, 



T-H 2 



(.B.43) 



Therefore, the local truncation error estimate for t± is obtained, by a combination 
of (|BT2|) , (|BT7|) , (IB~23| . (jB~24|) . (|B~30l) . (|B~43|> and a detailed comparison between 
the truncation equation (I3.22p , (|3.23l) and the original PDE: 



(38?$ + <9 t $ - A ($ 3 ) - (1 - e)A$ - 2A 2 $ - A 3 $ 

= PF\ en + F2 en — Fz en — (1 — e)F± en — Fc, en — i^en = 0. 



(B.44) 



In addition, the constant estimate (|3.25p for M is also satisfied, by a careful check. 
The estimate for t-i is very similar. We denote the following quantity 



F- 



k+l/2 _ 



= [d t <t> N + -dt$ N )(-,t 



ik+l/2\ 



(B.45) 



A detailed Taylor formula in time gives the following estimate: 



pfc+l/2 _ p/c+l/2 



fc+l/2 

r 21 , with 



T) 



(B.46) 



at each fixed grid point Meanwhile, it is clear that F7 has the following decom- 

position: 



F 



fc+1/2 



+ 



d t $ h N +1 + d t $ k N s 2 d 3 $^ +1 + d 3 $* 



12 



pfc+1/2 pfe+1/2 

r 7,l r 7,2 



(B.47) 
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To facilitate the analysis below, we define two more quantities: 

<2 1/2 — ^ t 3 <M ■ ,i fc+1/2 )- (B.48) 
A detailed Taylor formula in time gives the following estimate: 

pfc+l/2 pfc+1/2 _ fe+1/2 pfe+1/2 pfc+1/2 _ fe+1/2 . , 

*7,1 -*7e,l ~ T 22 . ^7,2 ~ ^7e,2 ~ T 23 > wltn 

l|r 2 2]] if(0)T) < Cs 3 \\*n\\hho,t) ^ Cs " ll*IU«(o,T) ■ (B.49) 
l|r23]Lj (0 ,T) < Cs 3 ||$iv|| ff 4 (0)T) < Cs 3 \\n H Ho,T) . (B.50) 

at each fixed grid point Consequently, a combination of (|B.46I) - (|B.50[) shows 

that 

F fe+i/2 _ ^fe+1/2 = T *H-V2 >with || T2 || i?(o t) < Cs 3 \m\ m(0<T) . (B.51) 
This in turn implies that 

11*2 - f 7\\li ( o,T:LI) < Cs 3 11*11^4(0,^) . (B-52) 

which is exactly (|3.23p . Also, the constant M associated with 72 satisfies (|3.25p . The 
consistency analysis is finished. 

REFERENCES 

[1] R. Backofen, A. Ratz, and A. Voigt. Nucleation and growth by a phase field crystal (PFC) 

model. Phil. Mag. Lett., 87:813, 2007. 
[2] A. Baskaran, Z. Hu, J.S. Lowengrub, C. Wang, S.M Wise, , and P. Zhou. Energy stable and 

efficient finite-difference nonlinear multigrid schemes for the modified phase field crystal 

equation. J. Cora-put. Phys., 2012. in review. 
[3] J. Boyd. Chebyshev and Fourier Spectram Methods. Dover, New York, second edition, 2001. 
[4] C. Canuto and A. Quarteroni. Approximation results for orthogonal polynomials in sobolev 

spaces. Math. Comp., 38:257-276, 1982. 
[5] W. Chen, S. Conde, C. Wang, X. Wang, and S.M. Wise. A linear energy stable scheme for a 

thin film model without slope selection. J. Sci. Comput., 2011. (accepted). 
[6] M. Cheng and J. A. Warren. An efficient algorithm for solving the phase field crystal model. J. 

Comput. Phys., 227:6241, 2008. 
[7] K.R. Elder and M. Grant. Modeling elastic and plastic deformations in nonequilibrium pro- 
cessing using phase filed crystal. Phys. Rev. E, 90:051605, 2004. 
[8] K.R. Elder, M. Katakowski, M. Haataja, and M. Grant. Modeling elasticity in crystal growth. 

Phys. Rev. Lett., 88:245701, 2002. 
[9] D. Eyre. Unconditionally gradient stable time marching the Cahn-Hilliard equation. In J. W. 

Bullard, R. Kalia, M. Stoneham, and L.Q. Chen, editors, Computational and Mathematical 

Models of Microstructural Evolution, volume 53, pages 1686-1712, Warrendale, PA, USA, 

1998. Materials Research Society. 
[10] Z. Hu, S. Wise, C. Wang, and J. Lowengrub. Stable and efficient finite-difference nonlinear- 

multigrid schemes for the phase-field crystal equation. J. Comput. Phys., 228:5323-5339, 

2009. 

[11] S. Gottlieb J. Hesthaven and D. Gottlieb. Spectral Methods for Time- dependent Problems. 

Cambridge University Press, Cambridge, UK, 2007. 
[12] V. Lebedev, A. Sysoeva, and P. Galenko. Unconditionally gradient-stable computational 

schemes in problems of fast phase transitions. Phys. Rev. E, 83, 2011. 
[13] J. Mellenthin, A. Karma, and M. Plapp. Phase-field crystal study of grain-boundary premelting. 

Phys. Rev. B, 78:184110, 2008. 



24 



A. BASKARAN, J.S. LOWENGRUB, C. WANG, AND S.M. WISE 



[14] J. Shen, C. Wang, X. Wang, and S.M. Wise. Second-order convex splitting schemes for gradient 
flows with ehrlich-schwoebel type energy: Application to thin film epitaxy. SI AM J. Numer. 
Anal, 50:105-125, 2012. 

[15] P. Stefanovic, M. Haataja, and N. Provatas. Phase-field crystals with elastic interactions. Phys. 

Rev. Lett., 96:225504, 2006. 
[16] Peter Stefanovic, Mikko Haataja, and Nikolas Provatas. Phase field crystal study of deformation 

and plasticity in nanocrystalline materials. Phys. Rev. E, 80:046107, 2009. 
[17] J. Swift and P.C. Hohenbcrg. Hydrodynamic fluctuations at the convective instability. Phys. 

Rev. A, 15:319, 1977. 

[18] C. Wang, X. Wang, and S.M. Wise. Unconditionally stable schemes for equations of thin film 

epitaxy. Discrete Cont. Dyn. Sys. Ser. A, 28:405-423, 2010. 
[19] C. Wang and S. Wise. An energy stable and convergent finite-difference scheme for the modified 

phase field crystal equation. SIAM J. Numer. Anal., 49:945-969, 2011. 
[20] C. Wang and S.M. Wise. Global smooth solutions of the modified phase field crystal equation. 

Methods Appl. Anal., 17:191-212, 2010. 
[21] S.M. Wise. Unconditionally stable finite difference, nonlinear multigrid simulation of the Cahn- 

Hilliard-Hclc-Shaw system of equations. J. Sci. Comput., 44:38-68, 2010. 
[22] S.M. Wise, C. Wang, and J.S. Lowengrub. An energy stable and convergent finite-difference 

scheme for the phase field crystal equation. SIAM J. Numer. Anal, 47:2269-2288, 2009. 



